Assessing wind field characteristics along the airport runway glide slope: an explainable boosting machine-assisted wind tunnel study

Aircraft landings are especially perilous when the wind is gusty near airport runways. For this reason, an aircraft may deviate from its glide slope, miss its approach, or even crash in the worst cases. In the study, we used the state-of-the-art glass-box model, the Explainable Boosting Machine (EBM), to estimate the variation in headwind speed and turbulence intensity along the airport runway glide slope and to interpret the various contributing factors. To begin, the wind field characteristics were examined by developing a scaled-down model of Hong Kong International Airport (HKIA) runway as well as and the surrounding buildings and complex terrain in the TJ-3 atmospheric boundary layer wind tunnel. The placement of probes along the glide slope of the model runway aided in the measurement of wind field characteristics at different locations in the presence and absence of surrounding buildings. Next, the experimental data was used to train the EBM model in conjunction with Bayesian optimization approach. The counterpart black box models (extreme gradient boosting, random forest, extra tree and adaptive boosting) as well as other glass box models (linear regression and decision tree) were compared with the outcomes of the EBM model. Based on the holdout testing data, the EBM model revealed superior performance for both variation in headwind speed and turbulence intensity in terms of mean absolute error, mean squared error, root mean squared error and R-square values. To further evaluate the impact of different factors on the wind field characteristics along the airport runway glide slope, the EBM model allows for a full interpretation of the contribution of individual and pairwise interactions of factors to the prediction results from both a global and a local perspective.

www.nature.com/scientificreports/ Constant wind speed is preferable from a takeoff and landing safety perspective. Consequently, wind speed is not a problem in and of itself. The issue arises, however, when the wind speed and/or direction fluctuate abruptly 2,3 . Wind shear is the term for this phenomenon. According to the International Civil Aviation Organization (ICAO), wind shear is defined as a sudden change in wind speed or direction within three nautical miles (~ 5500 m) of an airport runway threshold at an altitude of 500 m or less 4 . Changes in the aircraft's lift can be significantly impacted by wind shear. The declining and raising headwind shear on an aircraft during an approach is depicted in Fig. 2. When the approaching aircraft is subjected to a declining headwind, the lift is lowered, and the aircraft may possibly land short of the runway. In the case of a rising headwind, the aircraft's airspeed increases in relation to the surrounding air flow, generating more lift and resulting in a flatter angle of descent or even a climb, which ultimately may cause a missed approach 5 .
In addition to wind shear, turbulence may also affect an aircraft's landing. Rapid, irregular air motion is what causes turbulence. It causes jolts or hiccups. In extreme situations, the aircraft may lose control momentarily [6][7][8] . When an aircraft encounters turbulence, loose objects within the cabin may become dislodged, and passengers who are not wearing seat belts may sustain injuries 9 . A pilot must take corrective measures to ensure safety in the presence of significant wind shear and turbulence.
Causes of wind shear and turbulence. Two basic factors causing wind shear and turbulence are meteorological elements 10,11 and the geographical environment inside or around the airport [12][13][14] as depicted in Fig. 3. Tropical cyclones, thunderstorms, cold and warm fronts, and jet streams (narrow bands of strong winds) are common weather phenomena that can cause wind shear and turbulence. Sea breezes, strong monsoon winds, and strong winds that climb over hills and man-made structures are additionally recognized to cause low-level wind shear near runways. Similarly, the intensity of turbulence is expected to decline when the wind blows over flat terrain, such as open water, and to increase when the wind blows over sizeable obstacles, such as mountains and buildings 15,16 . Assessment of wind shear and turbulence by numerical simulation. Over   www.nature.com/scientificreports/ national Airport (HKIA), Lei et al. 17 used RANS equations and a LES model based on CFD to simulate terraininduced wind shear. Similar to this, Chen et al. 18 developed a high-resolution LES model using information from the Weather Research and Forecasting (WRF) model for the detection and analysis of terrain-induced low-level wind shear at HKIA. Boilley and Mahfouf 19 performed numerical simulations with the nonhydrostatic Meso-NH model to predict the low-level wind shear over Nice airport in France. Rasheed and Sørl 20 employed CFD analysis for the terrain-induced turbulence at Kristiansand airport, Kjevik. Zhang et al. 21 also developed CFDbased terrain-induced wind shear models for Beijing Capital International Airport (BCIA). Additionally, researchers have used CFD to assess turbulence intensity. Both RANS simulation and LES were used to study the transient nature of terrain-induced flow disturbances along the glide paths of descending aircraft 22 . Similar to this, Shimoyama et al. 23 used LES to gain insight into the turbulence near Japan's Shonai Airport. According to the aforementioned studies, it has been found that the CFD model has precisely modeled wind shear and turbulence in the close vicinity of the airport. However, these studies were naturally constrained in their duration and/or scope by the exclusive use of simulation models. They employed the RANS equation, which can simulate and forecast the average airport wind field characteristics but cannot gauge the actual wind shear and turbulence intensity.
Applications of wind tunnel studies. The wind tunnel studies offer a viable alternative to numerical simulation models for the evaluation of wind shear and turbulence characteristics surrounding airports. Consequently, the results of wind tunnel experiments can serve as a crucial basis for determining the reliability of numerical simulations. Numerous scholars from various fields have conducted wind tunnel experiments to evaluate the impact of wind on bridges, wind turbines, low-and high-rise building structures, etc. For instance, Diana et al. 24 investigated the wind-bridge interaction by employing wind tunnel experiments and considering the aerodynamic phenomena that affect the various bridge components, primarily the deck and towers. He and Zou 25 implemented a wind tunnel experiment to examine the aerodynamic problems that arise with train-bridge systems in the presence of wind. Li et al. 26 studied the wind characteristics at the bridge site in the deep-cutting gorge using wind tunnel experiments. In the simulated atmospheric boundary layer, the effects of varying incoming wind directions on the wind characteristics at the bridge site were investigated.
Similarly, Su et al. 27 demonstrated the influence of array configuration on the power performance of vertical-axis wind turbines through wind tunnel experiments. Tip speed ratio, rotational direction, and position were considered, and the performance of various array configurations was compared. Ice buildup on the blade surface of wind turbines installed in cold, humid climates significantly reduces their performance. Gao et al. 28 conducted icing wind tunnel tests to investigate the characteristics of icing under various tip speed ratios and rime ice conditions. Likewise, using wind tunnel experiments, Cai et al. 29 evaluated the impact of collective pitch control and cyclic pitch control on the power and loads of floating offshore wind turbines under the condition of uniform wind velocity.
Li et al. 30 analyzed the effect of winds on the wooden pagoda in China by employing wind tunnel tests and computing the wind loads on the pagoda. One of the parts of low-rise buildings that is most susceptible to damage in windstorms is the overhang on the roof. Therefore, Wang et al. 31 conducted field monitoring during typhoons and wind tunnel experiments to analyze the wind pressures on the upper and lower surfaces of the roof overhang in order to investigate the wind-induced pressure distributions and forces on the roof overhang of a typical low-rise building.
Despite the fact that a number of researchers have successfully employed wind tunnel experiments in their respective fields, the most serious drawback of wind tunnel investigations is the high cost of testing as well as the lack of testing infrastructure and time. To achieve the desired outcomes, numerous experiments are needed, carried out in many different kinds of settings. This reduces productivity and wastes both time and money. Empirical modeling techniques must replace experimental work in order to overcome the aforementioned limitations.
Combining machine learning with wind tunnel studies. In recent years, engineering has made great strides with machine learning strategies [32][33][34][35] . This is a result of the growing demand for sophisticated computational methods to process enormous data sets. Several researchers have utilized it to integrate machine learning strategies with the wind tunnel outputs, such as Weng and Paal 36 developed a novel machine learning-based www.nature.com/scientificreports/ wind pressure prediction model (ML-WPP) for non-isolated low-rise buildings by utilizing the non-isolated low-rise wind tunnel dataset. A machine learning-based predictive model of crosswind vibrations of rectangular cylinders was developed by Lin et al. 37 . Building pressure patterns were identified using an unsupervised machine learning method by Kim et al. 38 . Using the data collected from wind tunnel tests, deep learning methods to predict wind pressures on tall buildings were proposed 39 . Similarly, Wada et al. 40 used deep reinforcement learning in conjunction with discrete actions to control the pitch of an aerial vehicle based on experiments in a wind tunnel. Numerous applications have evidently led to the use of machine learning techniques to predict the windinduced responses of structures. In contrast, its applicability to addressing the impact of wind field characteristics on airport runway glide paths is severely constrained.
Research process. This research combines wind tunnel experiments with machine learning models to estimate wind shear and turbulence intensity along the airport runway glide path in order to overcome the shortcomings of numerical simulations. Nevertheless, one of the most significant drawbacks of machine learning models is their black box nature. They required an explicit post-hoc explanation tool for the interpretation of the model and factors. In this study, we propose an explainable boosting machine (EBM) strategy 41 combined with wind tunnel experiments to develop non-parametric models for estimating variation in headwind speed and turbulence intensity along the airport runway glide path. A scaled model of HKIA's north runway and Lantau Island (south of HKIA) was constructed in the TJ-3 atmospheric boundary layer (ABL) wind tunnel at Tongji University, Shanghai. Since HKIA is particularly susceptible to wind shear and turbulence due to its location and the topography of the area, an analysis of the airport's wind field is warranted 11,18 . The north runway at HKIA is oriented as 25RA and 07LA; consequently, the glide paths of 25RA and 07LA were subject to different inflow directions. Along the glide paths of both 25RA and 07LA, measurements were made at fixed distances from the runway threshold. The EBM model was trained and evaluated based on the results of wind tunnel tests to estimate variations in headwind speed and turbulence intensity. In addition, the performance of the EBM model in estimating variations in headwind speed, crosswind speed, and turbulence intensity was compared to the performance of other glass box models, including decision tree (DT) 42 and linear regression 43 , as well as its counterpart black box models, Random Forest (RF) 44 , Extreme Gradient Boosting (XGBoost) 45 , Adaptive Boosting (AdaBoost) 46 , and Extra Tree (ET) 47 . The Bayesian optimization approach 48 was employed during the training-validation stage to fine-tune the hyperparameters of machine learning models. Afterwards, EBM was also used to investigate a number of factors and their global and local pairwise interactions. Figure 4 depicts the entire process of wind tunnel experiments combined with machine learning modeling.

Data preparation
Wind field effect at HKIA. The HKIA is located on the island of Lantau, which is subtropical terrain, off the southeast coast of the Chinese mainland as shown in Fig. 5 49 . Numerous observational and modeling studies have shown that the HKIA's complex orography and high land-sea contrast are favorable for the emergence of low-level turbulence and wind shear [50][51][52] . Pilot flight reports based on the HKIA indicate that nearly 1 in 500 flights have been subjected to wind shear and turbulence since the airport's opening. 97% of the pilot reports revealed 15 to 25 knots of low-level wind shear. Terrain-induced wind shear was to blame for about 70% of the wind shear that pilots reported 53 . In addition to the terrain, nearby buildings are also significant sources of wind shear and turbulence 54 , as illustrated in Fig. 6 55 .
The wind tunnel experiments. This study employed wind tunnel experiments to evaluate the turbulence intensity and variations in headwind speeds along the glide path of the north runway at HKIA under various inflow wind conditions. At the State Key Laboratory for Disaster Reduction in Civil Engineering at Tongji University in Shanghai, experiments were carried out in the TJ-3 ABL wind tunnel. This wind tunnel is a low-speed, closed-circuit, vertical-return type. The wind tunnel's test chamber has dimensions of 14 m long, 15 m wide, and 2 m high. For the wind tunnel testing, Lantau Island, the north runway of the HKIA, and the buildings to the south of the runway, with a test range of 27.2 km and an average height of 425.2 m, were taken into consideration. The Lantau island and building were included in the 1:4000 scale model with a diameter of 6.8 m that was used in the wind tunnel for testing purposes, as shown in Fig. 7a. The height within the test range was 0.106 m. The surrounding terrain model was constructed layer by layer using dense foam with a one-inch texture (equivalent to a difference in actual terrain height of 40 m). The wind tunnel blockage ratio was calculated to be 2.402%, which is less than 5% (preferably for wind tunnel research) and satisfied the requirements of the wind tunnel experiments. The wind direction was varied from 90 to 240 degrees in 15 degree increments, taking into account the typical airflow patterns of east to southeasterly winds and the southwest monsoon of Hong Kong 56 , as shown in Fig. 7b. Therefore, eleven different wind conditions were in operation. The wind directions were given as follows: 0 degrees for the north wind, 90 degrees for the east wind, 180 degrees for the south wind, and 270 degrees for the west wind. Given that aircraft usually perform headwind landings, the intervention of a southwest wind (with wind directions of 190, 210, 225, and 240 degrees) was primarily taken into account for the glide slope of Runway 25RA, and the intervention of an east-southeast wind (with wind directions of 190, 210, 225, and 240 degrees) was taken into consideration for the glide slope of Runway 07LA.
The pink dotted line in Fig. 7b shows that aircraft typically follow a glide slope of 3 degrees in the final three nautical miles before landing on the runway 57 . In wind tunnel experiments, two sets of a total of eight measurement points (a1, a2, a3, a4) and (b1, b2, b3, b4) were placed along the glide paths of runways 07LA and 25RA, respectively. The cobra probes were mounted on bespoke stands, and the installation's overall height was adjusted so that it was in line with the height of the measurement site. The probes were oriented toward the inflow  www.nature.com/scientificreports/

Development of explainable boosting machine model
An explainable boosting machine (EBM) model was developed to estimate the turbulence intensity and variation in headwind along the airport runway glide slope using data from wind tunnel experiments. The EBM model is an interpretable model, in contrast to other machine learning models. When implementing machine learning models, interpretability can be crucial. By interpreting models, one can gain confidence in them and facilitate their adoption. It may also be useful for debugging the model, and in some cases, explanations for the model's predictions are required. Figure 8 depicts the EBM model, with factors serving as inputs and outputs. Prior to developing an EBM model, it is necessary to perform label encoding. The factors, which include the effect of buildings at or near HKIA, the orientation of the approach runway, and the distance from the runway, are coded as shown in Table 1.
Fundamental concept of EBM model. The EBM algorithm is based on GAMs, which are universally acknowledged to be the gold standard for comprehensibility. Given that = x j , y j is a training dataset with a length of R , that (x 1 , x 2 , ..., x R ) are the input vectors with v number of factors, and that y r is the target factor, then the GAMs takes the form as shown in Eq. (1).
where, x j indicates the jth factor in the factors set, is the link function that conforms the GAM to regression (e.g., = identity) or classification (e.g., = logistic), and Ŵ r is the attribute function. EBM, which employs bagging and gradient boosting, provides several significant advantages over conventional GAMs. An adequate description would be bagged, augmented, and lowered shallow trees. As depicted in the Fig. 9, the core of the algorithm employs boosted trees. In EBM, training is carried out repeatedly, with each repetition involving the creation of a distinct boosting procedure for each factor. The factor's order is irrelevant because of the low learning rates used in this process. In order to determine the best way to determine the contribution of each factor to the model's prediction, the high number of iterations aims to reduce the effects of co-linearity. Additionally, pairwise interaction terms    www.nature.com/scientificreports/ can be automatically detected and included by EBM, improving the model's accuracy while maintaining its explainability. The fact that EBM is an additive model allows for the capture and visualization of each factor's contribution, which improves explainability. Small trees are built sequentially for each iteration of the model, each of which can only use one factor. In a boosting manner, the residual (ε) is updated and a new tree is built using a different factor. Every iteration involves doing this for every factor. When training is finished, we can view all the trees that a particular factor used to build them. Due to the concepts of additivity and modularity, it is feasible to order and visualize the contributions to determine which factor had the greatest influence on the individual prediction. Another advantage of the EBM is the ability to improve accuracy by incorporating pairwise interactions into standard GAMs that becomes GA2Ms with the form and given as Eq. (2).
Here, the two-dimensional interaction Ŵ rs (x r , x s ) can be depicted as a heat map on the two-dimensional x r − x s plane while still retaining an exceptionally high degree of intelligibility. Initially, GA2M constructs the best GAM possible, and then it searches through the residuals to find all of the possible combinations of interactions, and then it ranks them based on how crucial they are. The last step is to integrate into the model the κ most significant pairwise interactions, where κ is a value determined through the process of cross-validation.
Hyperparameter tuning of EBM via Bayesian optimization. In this study, the process of fine-tuning the hyperparameters of EBM model as well as other counterpart black box and other glass box model was carried out with the aid of a Bayesian optimization strategy. It is one of the algorithms for global optimization, and it finds widespread application in the field of engineering 58,59 . The objective function that needs to be optimized in order to find the optimal values for the hyperparameters of the model can be expressed as Eq. (3).
where, x : Hyperparameters of the model, s : Hyperparameters search space, ϕ(x) : Objective function that reflects the association between models' performance and hyperparameters.
In this study, we employed R 2 as the evaluation metric. The purpose of this optimization is to identify the appropriate combination of hyperparameters x in such a way that the performance of the model ϕ(x) is increased to its full potential. The concepts of the Bayesian theorem are utilized in Bayesian optimization, which includes the following [Eq. If the Gaussian process is used, then the black box function ϕ will behave in a way that is consistent with the Gaussian distribution. At this very moment, the acquisition function ϒ(x) , which is based on the expected improvement, looks like this: where, : Cumulative distribution function of Gaussian distribution, β : PDF of Gaussian distribution The process of Bayesian optimization of hyperparameters, the iterative search is as follow: Iterate the loop ω times (ω = 1, 2, . . . , ) : 1. In accordance to the acquisition function ϒ(x) computed by Eq. (4), the hyperparameters combination x ω+1 = arg max ϒ(x; T) of the next cluster of instances is obtained. 2. Corresponds to an instance x ω+1 , compute the performance ϕ ω+1 of the machine learning model. Performance measures. To assess the performance of different models, four distinct metrics could be utilized, i.e., mean absolute error (MAE), mean squared error (MSE), root mean squared error (RMSE), and Coefficient of determination (R 2 ). The MAE is the average absolute value of the individual prediction errors across all instances [Eq. (7)]. MSE is calculated using the average squared difference between observed and predicted values, as shown in Eq. (8). RMSE is the square root of the difference between the observed and predicted www.nature.com/scientificreports/ values, according to Eq. (9). R 2 which ranges from 0 to 1, indicates the ability of a model to accurately predict values. R 2 is given by Eq. (10).
where, the total number of data points are denoted by ζ, the actual experimental outcomes are represented by σ x , the predicted outcomes are represented by σ x and σ avg represents the average of predicted outcomes.

Results and discussion
This section presents the development of the EBM model, other glass box models, such as DT and linear regression, and their counterpart black box models, such as RF, XGBoost, AdaBoost, and ET, to estimate variation in headwind speed and turbulence intensity based on data from wind tunnel experiments. In order to achieve better performance, it is crucial to note that the hyperparameters must be optimized prior to the development of these models. To achieve this, the results of the wind tunnel experiments were split into a training-validation set (comprising 70% of the total) and a testing set (comprising 30% of the total) 60 . The development of machine learning models and the tuning of their hyperparameters through Bayesian optimization were both done using the training-validation data. Bayesian optimization aimed to find the set of hyperparameters for various machine learning models that maximizes the R 2 value across a given sample space. Figures 10 and 11 show the progress of the Bayesian optimization process over 50 iterations using both glass box and black box models for variations in headwind speed and turbulence intensity, respectively. They were tuned using the hyperparameters associated with the models at the best R 2 to provide the best prediction model. Table 2 shows the optimal values for the hyperparameters that were obtained for the estimation of variation in headwind speeds, while Table 3 shows the optimal hyperparameters for models that estimate turbulence intensity.

Prediction of variation in headwind speed and turbulence intensity.
Once the optimal hyperparameters were determined, holdout or testing data was used to evaluate and compare the performance of various models. Importantly, the dataset was partitioned after being shuffled with 40% and 50% test data, and this demonstrated that the metrics values for model testing remained constant within the 95% confidence interval. In addition to tuning hyperparameters, this was performed to prevent over-fitting. When the size of the test data was altered, no anomalies in the performance metrics were observed. Table 4 presents the performance metrics of various models for estimating the variation in headwind speeds using both the training and test datasets. In terms of MAE (0.346), MSE (0.189), RMSE (0.435), and R2 (0.928) for the training dataset and MAE (0.571),   Similarly, in the case of both variations in headwind speed and turbulence intensity, we also reported and compared the prediction confidence and their ability to extrapolate on unseen data for both black box and glass box models using a residual scatter plot and a residual QQ (quantile-quantile) plot, as shown in Figs. 12 and 13. By analyzing the residual scatter plot's residual patterns, the model's applicability is shown. The discrepancy between the predicted and observed values is represented by residues. In a residual scatter plot, the target factor values are plotted on the horizontal axis, and the residuals are plotted on the vertical axis. When using the EBM model to analyze residuals for variations in headwind speed as well as turbulence intensity, in contrast to other models, the points are closer to the horizontal line at 0.00.

EBM interpretation.
Based on its performance metrics utilizing a testing dataset, EBM was inevitably the best model. From the EBM model, we can glean explanations. Consequently, we next utilized the EBM model to illustrate global and local interpretations of various individual and pairwise factor interactions.
Global importance of factors and pairwise interaction. The EBM algorithm is a generalized additivity model that enhances and expands the tree-based model. Due to the persistence of additivity, the contributions of factors can be ranked and plotted to demonstrate impact prediction in both global and local contexts. The exhaustive global explanation of the EBM makes it possible to visualize the effect of each possible combination of factors on the variation in headwind speed and turbulence intensity. Figure 14a and b provides a summary of the global significance of each individual factor and their pairwise interactions for variation in headwind speed and turbulence intensity, respectively. Wind direction contributed the most to the variation in headwind speed, followed by Runway. In pairwise interaction, the most significant contributor to the variation in headwind speed was the combination of Distance from Runway and Runway. Similarly, wind direction contributed the most to turbulence intensity, followed by runway length and distance from the runway. The combination of wind direction and runway contributed the most to the occurrence of intense turbulence in pairwise interaction. These results are also inline with the previous study 54,61 Additional output provided by the proposed EBM model for comprehending the global findings for the variation in headwind speed and turbulence intensity is shown as examples in Figs. 15 and 16, respectively. These Table 4. Performance measures of different models for estimating variation in headwind speeds.  www.nature.com/scientificreports/ include EBM shape functions and a heat map. For each predictor factor, the gradient boosting procedure yields a large number of decision trees, which are then used to generate the function characterizing the dependent factor's response to the specific predictor factor. One can then graphically represent these functions as a one-dimensional function, with the values of the predictor factors along the horizontal axis and the scores (i.e., the effect of the predictor factor on the prediction) along the vertical axis. When the score is greater than 0, it means that the value of the predictor factor is strongly linked to the outcome. In the case of the variation in headwind speed, Fig. 14a-c shows EBM shape functions and a heat map for the important individual and pairwise interaction factors. Figure 15a demonstrates that the variation in headwind speed is greatest when the wind direction is between 90° and 135°, as well as between 200° and 240°. The headwind speed along the glide path of Runway 25RA significantly decreases due to the presence of buildings at the airport that are in its path. Runway 25RA persists in being affected by a wind coming from the southwest at an angle ranging from 195° to 240°. Compared to Runway 07LA, Runway 25RA is more likely to experience a greater variation in headwind speed (Fig. 15b). This may be due to the presence of buildings close to Runway 25RA that cause variations in wind speed.
The pairwise interactions between runway and distance from runway are depicted in Fig. 15c. The heat map revealed that the variation in headwind speed is greatest within 0.5 miles of Runway 25RA. Consequently, shorter distances along Runway 25RA have been more susceptible to variations in headwind speed, and pilots may encounter severe wind shear in this area. Pilots must be extremely vigilant during the final approach to Runway 25RA. Figure 16a-c displays EBM shape functions and a heat map for the crucial individual and pairwise interaction factors in the case of the turbulence intensity. The higher turbulence intensity is more likely to occur for wind directions between 110° and 200°, as shown in Fig. 16a. Complex terrain in this direction blocks the wind, causing it to fluctuate in speed, which ultimately increases the intensity of the turbulence. In contrast to Runway 07LA, Runway 25RA is more likely to encounter higher turbulence intensity than Runway 07LA (Fig. 16b). This may be caused by the presence of building structures close to Runway 25RA that cause wind speed variations. This www.nature.com/scientificreports/ suggests that nearby structures are restricting the wind close to the runway 62 . The heat map in Fig. 16c shows how the runway and wind direction interact pairwise. It demonstrates that Runway 25RA is most likely to experience high turbulence intensity when the wind direction is between 170 and 220 degrees and that Runway 07LA is most likely to experience high turbulence intensity when the wind direction is between 130 and 170 degrees.
Local factor interpretation. The local explanation, which includes the explanations behind particular predictions in a single case, can be investigated for both variations in headwind speeds and turbulence intensity. The intercept is shown to be constant in gray, while additive terms that have a positive effect on output are shown in orange, and additive terms that have a negative effect on output are shown in blue, demonstrating local factor contribution assessments for each instance. The estimates for each instance are derived from the shape functions and pairwise interactions (heat map) based on the input values. We take into account two arbitrary cases where the predictions are very close to actual values, each for variations in headwind speed and turbulence intensity. Nevertheless, we can perform local interpretation in other instances as well. The local explanation for a variation in headwind speed prediction is shown in Fig. 17a. It demonstrates that the variation in mean headwind speed was high for a randomly chosen instance at Runway 25RA and a shorter distance from the runway threshold. Similarly, Fig. 17b shows the degree of turbulence for a randomly selected instance, showing that a wind direction of 240 degrees has a negative effect on the turbulence intensity. However, a closer proximity to the runway, the presence of buildings, and Runway 25RA all increase the likelihood of experiencing high turbulence.

Conclusion and future work
In order to estimate the variation in headwind speed and turbulence intensity along the airport runway glide slope based on wind tunnel experiments, this study employs an explainable boosting machine (EBM), a stateof-the-art glass box model. The proposed approach is both explicable and intuitive, and it provides accuracy on www.nature.com/scientificreports/ par with black box competitors. Conventional machine learning regression models lack explainability, but the EBM model has been effective in overcoming this limitation. Therefore, the following conclusions can be drawn: • Based on the testing data set for the variation in headwind speed, the Bayesian optimized EBM model had the best overall performance of all the black box and glass box models, with an revealed that "wind direction" and "Runway" significantly contributed to the occurrence of high variation in headwind speed. In the case of pairwise interaction, it was found that the combination of "distance from runway" and "runway" contributed more to the high variation in headwind speed. • The EBM-based interpretation also demonstrated that, in the case of individual factors, "Wind Direction" and "Runway" significantly influenced the occurrence of high turbulence intensity. It was found that the combination of "Wind Direction" and "Runway" contributed more to the high turbulence intensity in the case of pairwise interaction.
This proposed EBM framework could be applied to a comprehensive analysis of wind field characteristics at other airports as well. It is indeed a great resource for individuals associated with civil aviation. However, there are a few recommendations for future research.
• Although this study utilized a number of different input parameters to predict wind field characteristics along the airport runway glide slope, many other parameters, such as atmospheric pressure and temperature, could be considered in future studies. www.nature.com/scientificreports/ • In this study, the variation in headwind speed and turbulence intensity along the airport runway glide slopes were the parameters of concern. Similarly, the turbulence integral length scale as well as the variation in crosswind speeds at different points on the glide path can also be considered for future work.

Data availability
The datasets used and analyzed in the current study is available from the corresponding author on reasonable request.